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Abstract 

The Fermi-Pasta-Ulam (FPU) problem consists of the nonequiparti- 
tion of energy among normal modes of a weakly anharmonic atomic chain 
model. In the harmonic limit each normal mode corresponds to a peri- 
odic orbit in phase space and is characterized by its wave number q. We 
continue normal modes from the harmonic limit into the FPU parame- 
ter regime and obtain persistence of these periodic orbits, termed here 
q- Breathers (QB). They are characterized by time periodicity, exponen- 
tial localization in the g-space of normal modes and linear stability up 
to a size-dependent threshold amplitude. Trajectories computed in the 
original FPU setting are perturbations around these exact QB solutions. 
The QB concept is applicable to other nonlinear lattices as well. 

1 Introduction 

In the year 1955 E. Fermi, J. Pasta, and S. Ulam (FPU) published their cel- 
ebrated report on thcrmalization of arrays of particles connected by weakly 
nonlinear springs pQ. Instead of the expected equipartition of energy among 
the normal modes of the systems, FPU observed that energy, initially placed 
in a low frequency normal mode of the linear problem with a frequency io q and 
a corresponding wave number q, stayed almost completely locked within a few 
neighbour modes, instead of being distributed among all modes of the system. 
Moreover, recurrence of energy to the originally excited mode was observed. 
The door was thus opened to study the fundamental physical and mathemati- 
cal problem of energy equipartition and ergodicity in nonlinear systems, which 
involves the Kolmogorov-Arnold-Moser (KAM) theorem, dividing thresholds 
between regular and chaotic dynamics, and soliton bearing integrable models. 

From the present perspective, the FPU observation (equally coined FPU 
problem or FPU paradox) appears to consist of three major ingredients: (FPU- 

1) for suitable parameter ranges (energy, system size, nonlinearity strength) 
low frequency excitations are localized in g-space of the normal modes; (FPU- 

2) recurrence of energy to an initially excited low-frequency mode is observed; 
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(FPU-3) different thresholds upon tuning the parameters are observed - a weak 
stochasticity threshold (WST) which separates regular from chaotic dynamics, 
yet possibly preserving the localization character in q-space, and an equiparti- 
tion threshold (ET) also coined strong stochasticity threshold separating local- 
ized from delocalized dynamics in q-space (see E3 El for a review) . 

Two major approaches were developed. The first one, taken by N. Zabusky 
and M. Kruskal, was to analyze dynamics of the nonlinear chain in the contin- 
uum limit, which led to a pioneering observation of solitary waves 0. It took 
(FPU-1) as given, and aimed at obtaining quantitative estimates for (FPU-2). A 
second approach was proposed by F. Izrailev and B. Chirikov |5 who associated 
energy equipartition with dynamical chaos and aimed at an analytical estimate 
of the ET by computing the overlap of nonlinear resonances |7j which leads to 
strong dynamical chaos insuring energy equipartition. Below the threshold the 
dynamics is regular and, thus, no equipartition should occur. It aimed mainly 
at (FPU-3). Several other analytical 8, 9 and numerical 2D],5I]>ff2] threshold 
estimates have been published since, and will be discussed below. Note that 
similar effects have been observed in many other nonlinear discrete chain or 
field equations on a finite spatial domain, see e.g. |13j . 

In this work we show that stable periodic orbits, coined g-breathers (QB), 
persist in the nonlinear FPU chain, which are exponentially localized in q-space. 
The existence of these orbits was first reported in ref.^3]. Stability of these 
periodic orbits implies that nearby trajectories will be localized in q-space as 
well, and evolve in a nearly regular fashion for long times on submanifolds 
similar to low-dimensional tori in phase space. Recurrence times - to come 
close to an initial point on such a torus again - will be much shorter than the 
general recurrence times estimates for the dynamics on an assumed KAM torus 
in the whole phase space. That follows from the fact that a KAM torus in 
general will have the dimension equivalent to the number of degrees of freedom, 
and thus much larger than the effective dimension of the perturbed q-breather 
evolution. Upon tuning the parameters of the system, QBs will turn from stable 
to unstable at certain threshold values, allowing for low-dimensional chaotic 
evolution of nearby trajectories. The localization length in q-space depends on 
these parameters as well, and at critical parameter values QBs delocalize, leading 
to equipartition. The q-breather concept allows thus to address simultaneously 
all three FPU ingredients. At the same time it can be extended to completely 
different lattice systems and even to nonlinear field equations. 

In the present work we construct QBs continuing them from the linear case 
and study their properties both numerically and by an analytical asymptotic 
calculation. We compare the thresholds of QB existence and stability to the var- 
ious stochasticity thresholds mentioned above. Finally we show the persistence 
of QBs in thermal equilibrium and during long transient processes. 

2 The model 

The FPU system is a chain of N equal masses coupled by nonlinear springs with 
the equations of motion containing quadratic (the a-model) 

x n = (x n+ i - 2x n + x n -i) + a[(x n+ i - x n ) 2 - (x n - x„_i) 2 ] (1) 
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or cubic (the /3-model) 

X n = (X n+ l - 2x n + X n -i) + f}[(x n +l ~ X n ) 3 - (x n - Xn^if] (2) 

interaction terms, where x n is the displacement of the n-th particle from its 
original position, and fixed boundary conditions are taken xq = xn + i = 0. A 
canonical transformation 



q— 1 ^ ' 

takes into the reciprocal wavenumber space with N normal mode coordinates 
Q q (t). The equations of motion then read 

JV 

Qq+ UJ qQq = 1 = ^q^l^m^q.l^mQlQm (4) 

V 2 ( iV + 1 )(,m=l 

for the FPU-a chain and 



(3 

Qq +^qQq = ~ 2(jV + 1) ^ ^ q UJlUJ m UJ n C q .l 

,rn ,nQ lQ 7nQ n (5) 

for the FPU-/? chain 0. Here the coupling coefficients B q j tm and C q .i. myn are 
given by 

-B 9 ,i,m = ^(<5g±i±m,0 — 3q±l±m,2(N+l)), (6) 
± 

Cg,Z,m,n = ^ (<$q±Z±m±ri,0 — <5<j±i±m±n, 2(7^+1) — ^±i±m±n,-2(iV+l))- (7) 
± 

The sum in and 10 is taken over all 4 or 8 combinations of signs, respectively. 
The normal mode frequencies 



7T<7 

2(iV+l)' 



= 2sin otAi I Tv ( 8 ) 



are nondegenerate. 



3 A crosslink to discrete breathers 

The system of equations 0,© corresponds to a network of oscillators with 
different eigenfrequencies. These oscillators are interacting with each other via 
nonlinear interaction terms, yet being long-ranged in g-space. Let us discuss the 
relation of this problem to the wellknown field of discrete breathers ^3 ^3 E| • 
Neglecting the nonlinear terms in the equations of motion the q-oscillators 
get decoupled, each conserving its energy 

E q = \(Q 2 q +u 2 q Q 2 q ) (9) 
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in time. Especially, we may consider the excitation of only one of the g- 
oscillators, i.e. E q ^ for g = go only. Such excitations are trivial time-periodic 
and ^-localized solutions (QBs) for a = /3 = 0. 

This setting is similar to the case of discrete breathers (DB), that are time- 
periodic and spatially localized excitations e.g. on networks of interacting iden- 
tical anharmonic oscillators, which survive continuation from the trivial limit 
of zero coupling ^HJ. Notably, DBs exist also in FPU lattices ^7] and exis- 
tence proofs have been obtained as well [191 1201 |2"T] . The reason for the generic 
existence of DBs is two-fold: the nonlinearity of each oscillator allows to tune 
its excitation frequency out of resonance with other nonexcited oscillators. The 
case of a linear coupling on the lattice ensures a bound spectrum of small ampli- 
tude plane waves, and thus allows for the escape of resonances of a DB frequency 
and its higher harmonics with that spectrum. The spatial localization of DBs 
in such a case is typically exponential for short-range interactions [15| . Among 
the wealth of theoretical results we stress two here. First, if the coupling on the 
lattice is short-ranged but purely nonlinear, the DB localizes in space superex- 
ponentially [22] • Second, if the coupling is algebraically decaying on the lattice, 
DBs localize only algebraically as well 23 . 

Let us compare these findings to the g-breather setting. The nonidentical q- 
oscillators are harmonic, but with frequencies being nonresonant. That ensures 
continuation of a trivial QB orbit from zero nonlinearity into the domain of 
nonzero nonlinearity. Indeed, similar to the existence proof of discrete breathers 
by MacKay and Aubry [18j by continuation from the trivial uncoupled limit, 
we can start at zero nonlinearity and excite mode go to energy E qo , with all 
other modes at rest. By properly choosing the origin of time the solution will 
be time-periodic and time-reversible, i.e. Q q {t + T) = Q q (t), Q q (t) = Q q (—t). 
Following JH] we first observe that 

N 

U q (t) = t) q + LO q U q + — ===== ^ WqUlVmBq,l, m UlU m (10) 
V 2 ( N + 1 > l, m =l 

maps any choice of time-periodic time- reversible functions U q (t;T) onto time- 
periodic and time-reversible functions u q (t;T) where the period T is variable. 
For a = the particular choice U q (t;T) = Q q (t) and T = 2Tr/u qo ensures that 
u q (t;T) = 0, i.e. the map (|10l) of time-periodic and time-reversible functions 
(spanning a Banach space) contains a zero. In the next step we have to check 
whether the map Q10[l is invertible at a = 0. Following MacKay and Aubry [T%] 
invertibility is trivially obeyed for all q ^ qo in l[10() . provided m27r/T ^ u> q . 
This nonresonance condition is generically fulfilled for any finite lattice size 
N and a period T = 2it /ui qo + e where e is a small but finite number of the 
order 1/N. As for go itself, we conclude that invertibility is given provided we 
demand e.g. fixing the value U qo (t — 0;T) to a given nonzero number. With 
these possibilities we arrive at the full invertibility of (|10[l at a = 0. We may 
then conclude (and skip the full derivation which follows ^Sl) that the Implicit 
Function Theorem can be applied, and the zero of 110[) for a = will persist 
up to some nonzero value of a. Of course the same procedure can be applied to 
the FPU-/3 model. 

As for their degree of localization in g-space (if any) , the above results on DBs 
suggest that two mechanisms counteract: purely nonlinear interaction favours 
stronger than exponential localization, while long-range interaction tends to 
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delocalize the QB. We also recall that the seemingly simple problem of periodic 
motion of a classical particle in an anharmonic potential of the FPU type follows 
from the differential equation 

x = —x — ax 2 — (3x 3 . (11) 

Bounded motion at some energy E yields a solution which is periodic with some 
period T{E) =2ir/Q and can be represented by a Fourier series 

x(t) = j2 A « elknt ( 12 ) 

which leads to algebraic equations for the Fourier coefficients A k 

A k = k 2 n 2 A k - aJ^AkiAk-kx - P A kl A k2 A k _ kl _ k2 . (13) 

fcl fcl,fc2 

Note that the equations 1|13JI have similar properties as compared to ©,© - 
interaction between the Fourier coefficients is nonlinear but long ranged. Yet 
it is well known that the bounded solutions l|12fl to Ijllfl are analytic functions 
x(t) and thus the Fourier series coefficients A k converge exponentially fast with 

k [23. 

Let us summarize this chapter with stating that from a mathematical point 
of view g-breather existence and their properties can be treated in a similar way 
to the methodology of discrete breather theory. The two different types of ex- 
citations, localized in real space and localized in reciprocal g-space, have much 
in common when represented in their natural phase space basis choice. Both 
representations are connected to each other by a simple canonical transforma- 
tion, which is nothing but a rotation of the phase space basis. In the following 
we will use analytical and computational tools developed for discrete breathers 
and analyze the properties of g-breathers. 

4 q-Breathers in the a-FPU system 
4.1 Numerical results 

Let us consider the a-model and start with showing the evolution of the original 
FPU trajectory for a = 0.25, N — 32 and an energy E = 0.077 placed initially 
into the mode with go = 1. We plot the time dependence of the mode energies 
E q (t) in Fig^ for the first five modes. The period of the slowest (go = 1) 
harmonic mode is T\ — 1-njui\ w 66.02. We nicely observe slow processes of 
redistribution of mode energies, recurrences and also even slower modulations 
of recurrence amplitudes on time scales of the order of 10 5 . Note also that on 
time scales comparable to T\ all mode energies show small additional oscillations 
- and it is easy to see that they correspond to frequencies which are multiples 
of uj\ . The localization in g-space is also nicely observed, with the maximum of 
E$ being eight times smaller than that of E\ . 

In order to construct a g-breather, let us choose first a — 0, excite a normal 
mode with q = qo to the energy E qo = E and let all other g-oscillators be at rest. 
With that we arrive at a unique periodic orbit in the phase space of the FPU 
model. We expect that the orbit will stay localized in g-space at least up to some 
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critical nonzero value of a (and similarly for the /3-model 26 ) . We proceed with 
a series of successful numerical experiments continuing periodic orbits of the 
linear chain to nonzero nonlinearity. These orbits as well as their Floquet spectra 
can be calculated using well developed computational tools f° r exploring 
periodic orbits. We choose a Poincare section plane {x s — 0, x s > 0}, where s = 
[2(N + l)/qo] corresponds to an antinode of the mode Q qo . We map the plane y 
(all phase variables excluding x s ) onto itself integrating the equations of motion 
Q until the trajectory crosses the plane again: y n+1 — J : {y n ). A periodic 
orbit of the FPU chain corresponds to a fixed point of the generated map. 
As the initial guess we use the point corresponding to the go-th linear mode: 

x n (0) = ^tQ 9o (0) sin (fg) , x n (0) = 0. The vector function Q = f(y) - y 
is used to calculate the Newton matrix J\f = dQ(y)i/dyj. We use a Gauss 
method to solve the equations G{y) — Af(y — y ') for the new iteration y 1 and 
do final corrections to adjust the correct total energy E. The iteration procedure 
continues until the required accuracy e is obtained: — y||/||y|| < £ (we 

have varied e from 10~ 5 to 10 -8 ), where ||y|| = max{|j/i|}. We also note that all 
periodic orbits computed here are invariant under time reversal, which means 
that there exist times to when Q q (to + t) = Q q (to — t) for all q. Since mode 
velocities vanish at these times, a QB can be also computed e.g. by initially 
choosing all mode velocities to zero, and integrating until Q qo vanishes again at 
the same sign of Q qo as it was at the starting point. Then the fixed point of 
the corresponding map of the mode coordinates only suffices to obtain an exact 
periodic solution. 

We have used one of the original parameter sets of the FPU— a study a = 
0.25, E = 0.077, N = 32 [J to find stable g-localized QB solutions with most 
of the energy concentrated in the mode go = 1 (and added the cases go = 2, 3 
for comparison, see Fig l2l3fl . Note, that in (JHJ 2(N + 1) is always a multiple 
of 2, and in the particular case of N = 32 it is also a multiple of 3. Then, for 
qa = 2 and qo = 3 only modes with q being a multiple of go are excited. Let 
us discuss the properties of the found solutions in some detail. In contrast to 
the original FPU trajectories, the QB is characterized by mode energies being 
almost constant in time. The straight lines in Fig^ compare the mode energies 
on the QB with go = 1 with the FPU trajectory. The period of the QB solutions 
is very close to the corresponding period T qa = 2n/uj qo of the harmonic mode 
which is continued. During one QB period of oscillations a relatively small 
energy interchange between the modes (of the order of 2 percent) is observed 
(Fig|3Ja-c)). with frequencies which correspond to multiples of the QB frequency 
- just like the small fluctuations of the mode energies for the FPU trajectory 
mentioned above. There are well defined intervals during which the energy of 
nonlinear coupling energy E c — E to t — J2 q =i E q (Etot being the full energy of the 
chain) increases sharply yet being overall very small (between 1 and 3 percent 
of the total energy, see Fig0Jd)). In Fig|S]we show the time dependence of the 
first five mode coordinates on the QB with go = 1. The time reversal symmetry 
is nicely observed (although not used for construction) at times to ~ 17,49. 
Note also that while Q± oscillates predominantly with the main QB frequency 
VLqb which is close to wi, the second mode Q2 is dominated by 2£Iqb, the third 
one Q3 by 3£Iqb etc. 

We conclude the numerical results on the a-FPU case with QB solutions for 
go = 1, E = 0.077 and various values of a up to a = 0.8, which are shown in 
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FigEU We observe that the localization length increases with increasing a, so 
there is a clear tendency towards an cquipartition threshold. 

4.2 Estimating the localization length 

To demonstrate localization of a QB solution in q-space analytically, we expand 
the solution to J3J into an asymptotic series with respect to the small parameter 
g = a/y/2{N+l): 

OO 

Q 9 (t)=X>"QW(t). (14) 

n=0 

Inserting this expansion into we obtain equations for the variables Qq (i). 
The equation for the zero order n = reads 

Q(°)+^Q( o )=0, (15) 

and for n > 

N n-1 

Q { q n) +U, 2 q Q { q n) = -U, q "l"mB q lm £ Q| ,ll) Q^ 2) . (16) 

l,m— 1 ^1,2—0 

ni-\-ri2=n—l 

As the zero-order approximation we take a single-mode solution to (j!5[l with 
Q { q °\t)^0 only for 9 = g : 

Q^ 0) =<WA cos^ t. (17) 

Consider the first-order equations (case n = 1 in ljl6Jl ). The right-hand 
part in (|16fl contains only one nonzero term corresponding to I = m = qo, 
rii = ri2 = 0. The coefficient B q i m here is non-zero only for q = qo and q = 2qo 
(here we assume 2qo < N, so that the second Kronecker symbol in ijBJ always 
equals zero). Thus, in the first order the only variables different from zero are 

Q&\t) and Qg>(t). 

Similarly, for n = 2, provided 3go < -W, we get Q q (t) ^ for q = g , 2qo, 3go 
only. 

The above allows us to formulate the following proposition. 
In the (k — l)th order of asymptotic expansion, provided kqo < N, variables 
Q q k (t) differ from zero only for q = qo, 2qo, ■ ■ ■ , kqo'- 

Qf- 1 \t)=Q Vq^{q ,2q ,...,kq }. (18a) 

It means, that the first non-zero expansion term for a mode q = kqo is of the 
order k — 1: 

O£(*) = Vm<*-1. (18b) 

In Appendix ^ we prove this statement by the method of mathematical 
induction and approximate the first non-zero term for a mode q = kqo at kqo <C 
N: 

Qk^W = A ^o (co S ku qo t + O ((kqo/N) 2 )) , (19a) 
where ^ 

A kqo = -fe. (19b) 

UJ qa 
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Ignored expansion terms lead to shifting the QB orbit frequency and next- 
order corrections to its shape. 

Multiplying l|19|) by a k ~ 1 and inserting it into ©, we approximate the mode 
energies as 

E kqa = fcV-%o ' (1 + O ((kqo/N) 2 )) (1 + 0(a)) , (20a) 

where 

7- V (20b) 

n % 

Note, that although mode energies are not strictly conserved in time, their 
variation is small, being limited to the orders of magnitude indicated in l|20a|) 
in parentheses (cf. Fig0J). We arrive at an exponential decay in g-space dressed 
with a power law: 

? \ , ? ,„..,!„ ( E qo 



ln(£ 9 ) =21n(^) + ^ln 7 + ln =32. (21) 
9o 9o \ 7 / 

with exponent q^Xwy for the exponential part. The predicted decay 1)20(1 fits 
nicely with a computed QB with qo — 1 and a = 0.025 as shown in [14| . In 
Figs l2l3l we compare (|20f) with the numerical results for a — 0.25 and qo = 
1,2,3. While qi shows that the analytical results underestimate the degree of 
localization with increasing a, we note that it does even at these values a better 
for larger values of qo . 

The calculation described is expected to fail at 7 close to or greater than 
one, when l|20|l does not support g-space localization. At the same time the 
comparison with the numerical results shows that higher order corrections to 
the analytical decay law extend the region of QB localization. Thus we estimate 
a lower bound for the QB derealization threshold energy E loc as 

£' oc = 7r 4 a- 2 (iV + l)- 3 ^. (22) 

The scaling behaviour is in good agreement with an analytical estimation of the 
equipartition threshold in the a-model due to second order nonlinear resonance 
overlap, suggested in 0, which reads a 2 EN 3 « q^. 

Finally we stress an even closer relation between the localization of g-breathers 
within the used perturbation theory and the Fourier series convergence of an- 
alytic periodic functions ^ll [l .l)12 [l .l)l.'j |l . One arrives from the QB problem to 
these equations by simply assuming u> q = const. 



4.3 Stability of g-breathers in the a-FPU system 

An important problem is the stability of QBs. For computing linear stability 
of an orbit Q q (t), the phase space flow around it is linearized by making a 
replacement 

Q q = QS) + Zi (23) 

in the equations of motion (jSJ) and subsequent linearizing the resulting equations 
with respect to £ 9 . Orbit stability is then characterized by the eigenvalues of 
the Floquet matrix, which defines the linear transformation of small deviations 
£ g by the linearized equations over one period of the orbit. If all eigenvalues /ij 
have the absolute value one, the orbit is stable, otherwise it is unstable |15II16| . 



8 



All QB solutions presented in this section are linearly stable, i.e. all eigen- 
values of the Floquet matrix which characterizes the linearized phase space flow 
around the QB orbits reside on the unit circle. This is at variance with the case 
of the /3-FPU model, which will be discussed below. 



5 /3-FPU system 
5.1 Numerical results 

Launching an FPU trajectory by exciting a single low-frequency mode leads to 
similar observations as for the a-model. Again energy is localized in g-space on 
a few modes, sometimes coined natural packets |29| . which again persists for 
very long times. 

It is possible to construct QBs in the /3-FPU model, following the same way 
as described above for the a-model. 

For numerical calculation we use a modified scheme. The section plane 
is defined in the g-space as {Q qo — 0,Q qo > 0}. It is parametrized as r = 
{Qq,q 7^ <7o}- The even coupling potential and fixed boundary conditions en- 
able us to introduce an additional constraint Q q (t = 0) = 0, and the veloc- 
ity Q qo is obtained using the condition of energy conservation Q qo (t = 0) = 

■\J2E — J2 q ^i qo Qq(t — 0)- A s in the case of the a-model, the QB is searched as 

a fixpoint of the mapping f* l+1 = T(r"-). 

We obtain QBs which are exponentially localized in g-space (FigjZJ. The 
smaller /3, the faster is the decay of the energy distribution with increasing 
wave number q. Note, that due to the parity symmetry of the /3-model (Eq.J2J 
is invariant under x n — > — x n for all n) only odd g-modes are excited by the 
<7o = 3 mode and get coupled [2E]- This follows also from the coupling matrix 
0. 



5.2 Estimating the localization length 

In the analytical computation the solution to © is expanded in powers of a small 
parameter p = 2 (n+i) ■ ^ n ^ ne n "th order of expansion variables Qq differ from 
zero at q — qo, 3qo, . . . , (2n + l)qo only. Then the first non-zero expansion term 
for a mode (2n+l)go is of the order n. Using the same approach of mathematical 
induction as described in Appendix^ mode energies in a QB are approximated 
as follows: 

E { 2k+i )qo = X k E qo -(l + O (((2k + l)q a /N) 2 )) (1 + O(p)) , (24a) 

where 

9/? 2 £2 (N + l) 2 

*° ■ (24b) 

Again, time variation of mode energies is limited to the orders of magnitude 
indicated in l|24a|l . We arrive at a pure exponential decay 

ln(E q )=laE qo + h±-l)\nX (25) 
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with exponent (2qo) 1 In A. The predicted decay (|24|) fits nicely with the com- 
puted QBs in Figd 

The QB derealization threshold is then estimated as 

E qo = l^qlr\N + l)-\ (26) 

The scaling behaviour is in good agreement with an analytical estimation of the 
equipartition threshold in the /3-model due to second order nonlinear resonance 
overlap, obtained in which reads (3 EN w q$ . 

5.3 Stability of g-breathers for the /3-FPU model 

First, we study QB stability by computing the Floquet matrix numerically and 
diagonalizing it. In FiglHlwe plot the absolute values of the Floquet eigenvalues 
of the computed QBs versus f3 for different system sizes N. QBs are stable for 
sufficiently weak nonlinearities (all eigenvalues have absolute value 1). When f3 
exceeds a certain threshold two eigenvalues get absolute values larger than unity 
(and, correspondingly, another two get absolute values less than unity) and a 
QB becomes unstable. Remarkably, unstable QBs can be traced far beyond the 
stability threshold, and moreover, they retain their exponential localization in 
q-space (Fig[7J|. As (3 is increased further, new bifurcations of the same type are 
observed. 

To study QB stability analytically in the first-order approximation, we write 
down the QB solution in the form 

Q q (t) = S qqo Acostbt + O(p), (27) 

where uj is the QB frequency, slightly shifted from u qa due to nonlinearity. The 
residual term 0{p) includes corrections to the QB orbit shape by expansion 
terms of order one and higher. 

The first-order correction to the QB frequency is determined by secularity 
caused by resonant nonlinear self-forcing of the mode qq in the first order of 
expansion: 

Q%+<Q% = -*<Q?/ (28) 
This equation is identical to that of an isolated oscillator with cubic nonlinearity 
(Duffing oscillator). The well-known expression of nonlinear frequency shifting 
in the Duffing oscillator then yields 

^^{ l + ^T) +0i/) ) (29) 

Linearizing equations of motion Q around (|27|l according to l|23|) , we arrive 
at a Mathieu equation (see Appendix El, which finally leads to an estimation of 
the Floquet multipliers which leave the unit circle due to a primary parametric 
resonance and cause instability: 

\»nn I = 1 ± i7 ^^i?-l + 0(^) (30) 

where 

R= 6(3E qo {N+ l)/ir 2 . (31) 
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The bifurcation occurs at R = 1 + 0(1/N 2 ). This instability threshold coincides 
with the criterion of transition to weak chaos reported by De Luca et al. [Sj. 
Note, that iffil does not contain the principal mode number g . Below the 
stability threshold (except for possible small high-order resonance zones) a set 
of stable QB modes exists. In the thermodynamical limit N — > oo, however, 
the energy of stable QBs tends to zero. Note that according to our analytical 
( Appendix lB|l and numerical results the instability modes correspond to q' — 
go i 1 and are even modes if go is odd and vice versa. The observed instability 
is thus connected to a lowering of the symmetry of the higher symmetry QB. 
That has been also observed to be the driving pathway for the onset of low- 
dimensional stochasticity in the FPU trajectory at the weak chaos transition 
|Sl,|nni> where the FPU trajectory acquires chaotic components in the time 
evolution, while still being localized in g-space. 

The result JHDJ is plotted in FigED with solid lines for N = 32, 64 and 
128, demonstrating good agreement with the numerical results. The agreement 
improves with increasing N |31| . 

Driscoll and O'Neil studied the instability of a single soliton in the continuum 
mKdV limit of the /3-FPU model |32| with periodic boundary conditions. A 
stability threshold obtained within the mKdV equation will qualitatively or 
semiquantitatively agree with the correct value obtained for the discrete chain, 
if the instability sets in for QBs which do not contain significant short wavelength 
components. The relation between single soliton and QB stability is less clear, 
since in the limit of vanishing nonlinearity QBs in chains or field equations with 
periodic boundary conditions correspond to standing waves, while single solitons 
transform into plane (running) waves. 

We conclude this section with the observation that the stability l|31|) and 
derealization l|26(l threshold estimates for QBs contain a unique parameter 
(3EN where E is the total energy, at variance with scaling estimates for the 
transition times to equipartition |33j . which obtain f3E/N instead. 

6 QBs and the FPU-trajectory for the a-FPU 
model 

Departing from the the QB orbit in phase space in the direction of the initial 
condition of the FPU trajectory implies adding to the QB solution, which is 
localized in g-space, a perturbation, which is localized in q-space as well. The 
perturbed trajectory will evolve essentially on a low-dimensional torus in phase 
space, whose dimension will correspond to the number of modes excited on the 
QB orbit - e.g. for a = 0.25 and go = 1 about four or five. In Fig|U]we compare 
snapshots of displacements at different times obtained for the original FPU tra- 
jectory in £Q and for the numerically exact QB solution from FigElfor a = 0.25 
and observe similar evolution patterns (see also Ref. 52)- Moreover, we took 
a series of points on a line which connected initial conditions of the FPU tra- 
jectory (E q9 Li — 0) with the numerically exact QB solution from Fig|2 For 
each of these points we integrated the corresponding trajectory and measured 
the average deviation A from the QB orbit. The dependence of A on the line 
parameter turns out to be an almost linear one, starting from zero when being 
very close to the QB orbit, and ending with a maximum value when being close 
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to the FPU trajectory. That supports the expectation that the FPU trajectory 
is a perturbation of the QB orbit. The FPU recurrence is gradually appearing 
with increasing A and is thus directly related to the regular motion of a slightly 
perturbed QB periodic orbit, which we tested also numerically. In Fig [Till we 
plot the mode energies E\ and E2 during integration of the FPU trajectory for 
a = 0.25 and qo = 1 after consecutive periods of the g-breather solution. The 
regular pattern indicates regular motion, and the thick dot, which corresponds 
to the g-breather solution itself, resides inside the quasiperiodic cloud of the 
FPU trajectory, indicating once more that the FPU trajectory is a perturbation 
of the q-breather and evolves around the QB in phase space. Zooming the time 
dependence of the mode energies for the FPU trajectories on time scales com- 
parable to the QB period shows very similar nearly periodic fluctuations as in 
FigQ] which are generated by the QB period. Using the linearized phase space 
flow around a QB we can estimate an effective recurrence time, which for the 
original FPU case is two times smaller than the recurrence time for the FPU 
trajectory. We tracked the change of the recurrence time with increasing A. 
When coming closer to the FPU trajectory, simply every second recurrence as 
observed for small deviations from the QB is suppressed, leaving us with the 
FPU recurrence time. That effect may be due to additional nonlinear contribu- 
tions to the phase space flow around a QB. 

Finally we note that extremely long computations of the FPU trajectory have 
been reported recently [23| - The trajectory localizes in g-space (and thus stays 
close to a g-breather) for times up to 10 10 . Only after that a mixing of mode 
energies is observed, possibly due to Arnold diffusion. Notably that critical 
time has been estimated by numerical scaling analysis for shorter transition 
times |12| . 

7 Towards transient processes and thermal equi- 
librium 

Once the existence and stability of QBs as exact solutions are established, it is 
interesting to analyze the contribution of these trajectories to the dynamics of 
transient processes and thermal equilibrium. It is well-known, that in states, 
corresponding to energy equipartition between degrees of freedom (or thermal 
equilibrium) energy distribution demonstrates strong statistical properties, with 
no energy concentrations on average in some subparts of phase space. This cir- 
cumstance, however, does not forbid existence of finite-time energy localizations 
whose lifetime may substantially exceed the characteristic period of plane waves. 
We recall the studies of discrete breather contributions to the dynamics of non- 
linear lattices in thermal equilibrium and transient processes |34j . 

Here we present results of numerical simulations of the /3-FPU chain (N = 
100, fixed boundary conditions) with two types of initial conditions: (i) all 
energy is located in a single mode E to t — E 3 — 1.58, E q ^ 3 = 0, (ii) all energy 
is randomly distributed among all modes: Q q (0) = £,k/^ q , Q q (0) — r/k, where 
£fc,77/c are random numbers, uniformly distributed in [— c, c], c taken to ensure 
E tot /N = 0.2. 

In the first case we take three values of f} = 0.6, 1.25, 5.0, for which energy 
delocalizes and essentially redistributes among the degrees of freedom after some 
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transition time T tr fFig lll|) . We note that with the chosen values of qo,Eo,N 
our previous results suggest that the QB becomes unstable at p » 0.01 and 
delocalizes at p « 1.6. The figures show a time window width 10 5 which is 
much larger than the largest QB periods which are of the order of 400. 

We observe that stochastic motion and energy flow to higher modes lead the 
system into a possibly long transient regime, in which QB-like objects (finite- 
lifetime single-mode excitations) are observed for all q (Fig lllf a.fr) , note that at 
different T tr qualitatively the same picture is observed). These objects survive 
for up to a hundred periods of phonon band oscillations. Since the smallest 
chosen values for P exceed the QB instability threshold value, we conclude 
that the QB instability is of local character in (/-space and does not carry a 
perturbed trajectory far away. When p is increased, the lifetime of higher 
frequency QBs drops down and they disappear in the high- and middle-frequency 
regions fFig lllf c). p — 1.25), and then remain observable only in the few lowest 
modes (Tig lllf d). P = 5.0). Note that for these values of p we already exceed 
the derealization threshold estimate for QBs. The observation of surviving 
low-g QB-like structures suggests that the energy flow between low and high q 
modes is sufficiently weak, so that some excess of energy is transferred into the 
large q domain, allowing for long-time energy localization in the low q domain. 

Similar effects can be observed for the second type of initial conditions 
(FigC2J) which mimics thermal equilibrium. At low P QBs can be observed in 
the whole phonon band frequency domain fFig ll2f a.fr. P — 0.05). As P is in- 
creased, high- and medium-frequency QBs disappear (Fig^fc,d), P — 0.1, 0.4). 
Note, that by rescaling coordinates and momenta (and, hence, energies) the 
set of total energies and nonlinear coupling strengths for both types of initial 
conditions practically coincide, so one may directly compare all subfigures from 
FigEUand FiglTJ 

8 Conclusion 

We report on the existence of g-breathers as exact time-periodic low-frequency 
solutions in the nonlinear FPU system. These solutions are exponentially local- 
ized in the g-space of the normal modes and preserve stability for small enough 
nonlinearity. They continue from their trivial counterparts for zero nonlinear- 
ity at finite energy. The stability threshold of QB solutions coincides with 
the weak chaos threshold in [H]. The derealization threshold estimate of QBs 
shows identical scaling properties as the estimate of equipartition from second 
order nonlinear resonance overlap UJ. Persistence of exact stable QB modes is 
shown to be related to the FPU observation. The FPU trajectories computed 
50 years ago are perturbations of the exact QB orbits. Remarkably, localiza- 
tion in g-space persists even for parameters when the QBs turn unstable. The 
QB concept - tracing time-periodic and q-space localized orbits together with 
their stability and degree of localization - allows to explain quantitatively and 
semiquantitatively many aspects of the 50 year old FPU problem. Moreover we 
show that dynamical localization in g-space persists for transient processes and 
thermal equilibrium remarkably well. Note that for certain cases specific sym- 
metries allow to obtain g-breathers with compact localization [2E1 OH! > where 
they have been interpreted as anharmonic plane wave solutions. 

The concept of QBs and their impact on the evolution of excitations in the 
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FPU system is expected to apply far beyond the stability threshold of the QB 
solutions reported in the present work. Generalizations to higher dimensional 
lattices and other Hamiltonians are straightforward, due to the weak constraint 
imposed by the nonresonance condition needed for continuation. QBs can be 
also expected to contribute to peculiar dynamical features of nonlinear lattices 
in thermal equilibrium, e.g. the anomalous heat conductivity in FPU lattices 

m 

A quantization of QBs will lead straightforwardly to quantum QBs, which 
will not differ much from their classical counterparts, due to the absence of 
discrete symmetries which may lead to tunneling effects. Thus a quantum QB 
state will be localized in g-space pretty much as a classical QB, allowing for 
ballistic-type excitations in the low-q domain. 

Finally we remark that the QB concept is not even constrained to an un- 
derlying lattice in real space. What we need in order to construct a QB is a 
discrete spectrum of mode energies, and nonlinearity which induces mode-mode 
interactions. Thus QBs can be constructed in various nonlinear field equations 
on a finite spatial domain as well, since one again obtains in the limit of small 
amplitudes a discrete mode spectrum, with the only difference to a finite lattice 
that the number of modes is infinite. 
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A QB localization in the a-FPU model 

First, we note, that (|18|l and l|19(l are true for k = 1, according to l|17|l . Further, 
we assume it is true for k = 1, . . . , n — 1 and prove, that it is then true for k = n 
as well. 

Equation l|16(l for the order n — 1 can be written as 

n-l N 

Q( n ~V + ^Q(«-D =-cj q J2 E ^mB^Q^Qt 2 ^ (32) 

ni,2=l l,m—l 
ni+n 2 =n 

Consider, for which mode numbers q the right-hand part in 1)32(1 is non-zero. As 
^1,2 < n, we obtain from (|18a|l Qi" 1 ^ for I = qq, 2qo, • • ■ , «i9o only, and 
Q^r 1] ? for 

m — <Z0i 2(^0 j ■ ■ ■ , n 2Qo only- According to the definition of B q i rn 
(JHJ) and taking into account that nqo < N, we obtain B q i m ^ for q = I + m 
and q = \l — m\. Then, the maximal mode number q, for which the right-hand 
part in 1-S2I) is non-zero, equals q ma x = "-i9o + ^2^0 = n Qo- This is equivalent to 
(I18a|l at k = n, which concludes the induction for (JTSJ. 
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To prove l|19|l . we write down l|32l) for q = nq^: 

n-l 

+ <oQi"* ^ = -Wn* E ^i^OfcW ^ ( 33a ) 

"1,2 = 1 

ni+n 2 =n 

Here we have taken into account, that the only non-zero term of the sum over 
I and m corresponds to I = n\qo, m — n,2qo, and that B nqQtniq0tn2qo = 1. 

As Ui,a < n , the variables Qi" 1 1 and q!™ 2 1 ' ) are expressed via Ijl9a(l . Then 
(I33a|l becomes a forced harmonic oscillator equation with a set of cosinusoidal 
harmonics in the right-hand part: 

n 

Ql^ + <oQ { n^ l) = E C ™ cosm^ i (33b) 

where C m are amplitudes of the harmonics. 

Omitting the free-motion part of the solution, we write down the expression 
for the forced oscillations of Qn qo ^ : 

n „ 

Qmo^ = £ u 2 ^ (34) 

m =0 "90 W qo 

Note, that the resonance denominator in (|34(l is minimal at m = n. In other 
words, driving frequency noj qa is the closest to the resonance among all the 
harmonics. Expanding the sine function in ||5J| into a Taylor series as sin a; = 
x — x 3 /6 + 0(x b ), we approximate the mode frequencies as 

c g = 4(V) 2 -i(g X ) 4 + 0((gx) 6 )) where * = ^ -> (35) 

The mentioned minimal denominator is then expressed at nqo/N — > as follows: 



2 2 2 

wl„ — n oj„ = 

nqo qo 22 



L»V- V.+O ((£)*) (3«) 



At the same time the relation of this denominator to any other one with m n 
is estimated as 



< -mM ° l^J ' (37) 

It means, that the harmonic m = n dominates in the solution l|34|l . while all 
the other harmonics are small values of the order O [(nqo/N) 2 ) with respect to 
this dominating one. 

Thus, (|34|l can be rewritten as 

«fe 1> -35^- e -^ , ( 1+0 (^)')) (38) 

Inserting l|19a|l (valid as ni,2 < n) into l|33a|l . we obtain the following ex- 
pression for the nth harmonic amplitude C n in i|33bjl : 

C n = -w ngo E w niW„ 2 ~A niqo A n2qo ■ (l + O ((Jjf) J J (39) 

"1,2 = 1 ~ V V / / 

rii+«2=« 
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Inserting l|39|) into (|38|) and taking into account (|35|1 . we arrive at 

Qn q o 1] = A nqo (cosruj go t + O ((^) 2 )) (40) 



where 

6 " 1 

An 1» = 2( n 2 _ 1 W, ^2 n l n 2 A niq A n2 q (41) 

^ > 10 m, 2 =l 

Expressing A ni9o and A„ 29o via l|19b(l (valid because m,2 < n) and taking 
into account that 

n— 1 _^ 

«i«2 = g« 2 (" 2 - 1), 

"1,2 = 1 

we obtain 

A™ 

A* = -^r- (42) 

Equations (|40|l and (|42(l coincide with (|19h ,.b^l at k = ii, thus concluding the 
induction. 



B QB Stability in the /?-FPU model 

Making a replacement l|23|) in the equations of motion JSJ and keeping only 
linear in £ g terms we obtain an equation describing the dynamics of infinitesimal 
deviations from the QB orbit: 

N 

£ q + oj% = -3pQ 2 go (t)u>l + 0( P 2 , 6), (43) 

r=l 

where 6 gr = uj q uj r C q0yq0yqyr , and 0{p 2 , £;) denotes a linear form of {£/} with small 
coefficients of the order 0(p 2 ). 

Inserting here Q qo (t) from J37J, we obtain the following Mathieu equation: 

N 

4 + U% = ~h(l + COsffi) J2 \r£,r + 0(h 2 ,^), (44) 

r=l 

where h = 3pE qol il = 2uj. 

In the vector-matrix form it can be rewritten as follows: 

£ + A£ + h(l + cosnt)B£ = 0(h 2 Z), (45) 

where £ = (£ g ) is a vector, A = (a qr ) is a diagonal matrix with elements 
2 B = (b qr ) is the coupling matrix. 

Further, we analyze parametric resonance in the equation (|45|l . treating h 
and Q as independent parameters, and then recall their dependence. 

In the limit h — > the equilibrium point £ = is strongly stable for all 
values of fl except for a finite number of values Sl n ki which satisfy 

uj k + loi = nVt nk i (46) 
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where n is a natural number, and the modes k and I belong to the same con- 
nected component of the coupling graph whose connectivity is defined by the 
matrix B. Strong stability implies that this point is also stable for all hamilto- 
nian systems close enough to the considered one. 

Each point (£l n ki,0) on the plane of parameters (Q,h) is associated with 
a zone of parametric resonance. Restricting ourselves to the case of primary 
resonance, which necessarily requires n = 1 in l|46|) . we specify the frequency £1 
as 

(l = (u k +ui)(l + S), (47) 

where the detuning parameter 6 is assumed to be the order 5 = 0(h). We seek 
for a solution to (|45|l in the form 

£= f m e ( - iQl+z+imQ)t + c.c. (48) 

m— — oo 

where ui\ = uif.(l + 5), f m = (/™) are unknown complex vector amplitudes, 
and z is a small unknown complex number. We make an assumption z = 0(h) 
which is confirmed further in the course of the calculation. 

Inserting (|48|l into (|45|l we obtain a system of algebraic equations for the 
amplitudes /™: 

[2iz(Qx + mil) - (wi + mO) 2 + w*] /"+ 

+ (hBf m ) q + t^B(f m+1 + f^ 1 )) = 0(h 2 \$\) (49) 

Note, that if the coefficient in square brackets in (|49|) is not a small value 
of the order 0(h), then the corresponding amplitude /™ is itself a small value 
of the order 0(/i|£|). Let us find out, for which values of the indices m and 
q the mentioned coefficient is small. As we assume z = 0(h), the difference 
— (uji + mil) 2 + uj 2 q needs to be a small value. It will be so, if the absolute 
value |d>i + mf2| is close to one of the eigenfrequencies uj q . According to the 
definition of Cj\ and the expression for £1 (|47|) . this implies \(m + l)i>Jk + rncoi \ = 
uj q . Gcnerically, due to the incommensurate eigenfrequencies spectrum, this 
condition is only fulfilled for m = 0, q = k or m = —1, q = I. 

It follows, that all amplitudes f™ except for f% and are small values of 
the order 0(h\£\). Then we are able to write down a closed system for /° and 
f^ 1 accurate to 0(h 2 \£\): 

(2izu> 1 + 2u; 2 (h ~ 8))f k + 1 = 0(h 2 \Z\) (50a) 

^b lk f k + (-2izCj 2 + 2uj 2 (h - 5))f^ = 0(h 2 \£\) (50b) 

where LJ2 = — = W;(l + 5). Note, that for a primary resonance the coupling 
coefficient bki must be non-zero, which means that the mode oscillators k and I 
are directly coupled. 

A nontrivial solution to this system exists if the determinant of the left-hand 
part (with an error 0(h 2 ) allowed for in each element) equals zero. From this 
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condition we derive 



(51) 

In the initial problem both h and f2 depend on the nonlinearity magnitude 
(3E qo . This dependence defines a line starting from the point (2w mo ,0) on the 
(fi, h) plane. The intersections of this line with the resonance zones are the 
regions of the QB orbit instability. 

The nearest primary resonance corresponds to k = too — 1, I = too + 1, 
n = 1. In this case we can derive a simpler expression for in the vicinity of 
the bifurcation point (near the edge of the resonance zone) if we let N — > oo at 
the same time assuming q$/ (N + 1) = const which means uj qo — const. 

From l|35(l we obtain 

uj qo+1 - uj qo -i = 2coaq >c- (1 + 0(k 2 )) (52a) 
Uto+xUn-i = < • (1 + 0{x 2 )) (52b) 
^o+i + = 2w go • (1 + 0{y?)) (52c) 

Then we express S from l|47|l : 

5= \ = -h+\>? + 0{x l ) (53) 

u qo +i + u qo -i 4 2 

Inserting l|52h -c) and H53JI into (|51|l and taking into account that if we are 
close to the bifurcation point (the expression under the square root in (|51(l is 
close to zero) then h — 0(x 2 ), we find 

1 

: 2 



zi, 2 = ±-w qo »?y/R - 1 + 0(x 2 ) + i ■ O(^) (54) 



where R = h/x 2 = &f3E qo (N + 1) /vr 2 . 

According to (|48|l . the absolute values of the Floquet multipliers involved in 
the resonance are calculated as 



27rRezi, 2 , , . 

\*hh I = ex P( q h ( 55 ) 



which finally leads to (|30() . 
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Figure 1: Evolution of the linear mode energies for the first five modes on a large 
timescale for (i) the original FPU trajectory for a = 0.25, E = 0.077, N = 32 
PQ (oscillating curves) and (ii) the exact QB solution (almost straight lines). 
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Figure 2: Snapshots of the linear mode coordinates Q q (symbols) along with an- 
alytical predictions (lines) for QBs with q = 1, 2, 3 for a = 0.25, E = 0.077, N = 
32, at the moment when all velocities Q q equal zero. Some modes have zero 
contributions and their corresponding symbols are not plotted here. 
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Figure 3: Distributions of the linear mode energies E q in q— space for QBs with 
qo = 1,2,3 for a = 0.25, E = 0.077, N = 32. 
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Figure 4: Evolution of the linear mode energies (a) E\, (b) E2, (c) E3 and the 
energy of nonlinear coupling (d) E c (see the text for definition) for the QB with 
q = 1, a = 0.25, E = 0.077, N = 32. 
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Figure 5: Evolution of the linear mode coordinates (51,2,3,4,5 for the QB with 
q Q = 1, a = 0.25, E = 0.077, N = 32. 



23 



LU 




5 10 15 20 25 30 



q 

Figure 6: Distributions of the linear mode energies E q for a = 0.1, 0.25, 0.4 and 
0.8 for QBs with q Q = l,E = 0.077, N = 32. 
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Figure 7: Energy distributions between g-modes in QBs for different nonlinear 
coupling coefficients versus q in linear and log scales with analytical esti- 
mations of the QBs exponential localization (dashed lines). Parameters are 
E = 1.58, go = 3, N = 32. Only odd modes are shown (see text). The symbols 
for g =/= 3, 9, 15, 21, 27 represent upper bounds, the real mode energies might be 
even less. Note that QBs persist even far beyond the stability threshold (see 
FigEJl. 
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Figure 8: Absolute values of Floquet multipliers of QBs with the energy 
E = 1.58 and qo — 3 and different N versus f3. Symbols: numerical results, 
lines: analytical results. 
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Figure 9: Snapshots of displacements (a) of the original FPU trajectory for 
a = 0.25, E = 0.077, N = 32 pQ and (b) of the corresponding exact QB solution 
from FigI3 taken at different times. 
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Figure 10: Map of the mode energies E\ and i?2 on the FPU trajectory after 
consecutive periods of the corresponding q-brcather for N = 32, a — 0.25 and 
<7o = 1. The thick dot is the results for the QB orbit. 



27 



x 10 



1.06 



1.04 



(a) 

i Vi. 



'if , : ; i. ''. s b J 



n 



20 40 q60 80 



X10 7 (C) 



1.008 

,.y :■ ••• 
If : | 



1.006 r ■' 

1.004 ' . v. ' . ' S 

; i #|| ! >. . . - v/. - 

1.002 



0.18 

0.16 

0.14 

0.12 

0.1 

0.08 

0.06 

0.04 

0.02 



0.18 

0.16 

0.14 

0.12 

0.1 

0.08 

0.06 

0.04 

0.02 



20 40 q 6c HO 



x 10 

I 13 



(b) 



1.006 

I . . i : :. •. 1 



1.004 
1.002 
1 



1.008 ,- 



20 40 q 60 80 100 

(d) 



1.006 .• ' 

1.004 ; [•_■■ 



1.00? .' 



20 40 q 60 80 100 



10.15 



0.05 



0.18 

0.16 

0.14 

0.12 

0.1 

0.08 

0.06 

0.04 

0.02 



Figure 1 1 : Space-time plots of mode energies E q evolving from the initial local- 
ization in the third mode for N = 100, E 3 {0) = 1.58 and (a) (3 = 0.6, T tr = 10 6 , 
(b) (3 = 0.6, T tr = 10 7 , (c) (3 = 1.25, T tr = 10 7 , (d) (3 = 5.0, T tr = 10 7 . 
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Figure 12: Space-time plots of mode energies E q evolving from random initial 
conditions for N = 100, E tot /N = 0.2 and (a) (3 = 0.05, T tr = 10 5 , (b) (3 = 
0.05, T tr = 10 6 , (c) = 0.1, T tr = 2- 10 5 , (d) j3 = 0.4, T tr = 2 ■ 10 5 . 
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